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A SCHEME FOR COMPUTING SURFACE LAYER TURBULENT FLUXES 
FROM MEAN FLOW "SURFACE OBSERVATIONS" 

MARTIN I. HOFFERT* 

Department of Applied Science , New York University , New York, NY 10003 , USA 

and 

JOEL ST0RCH + 

Brookhaven National Laboratory , Upton, NY 11973, USA 

Abstract. A physical model and computational scheme are developed for 

\ 

generating turbulent surface stress, sensible heat flux and humidity flux 
from mean velocity, temperature and humidity at some fixed height in the 
atmospheric surface layer, where conditions at this reference level are 
presumed known from observations or the evolving state of a numerical atmo- 
spheric circulation model. The method is based on coupling the Monin Obukov 
surface layer similarity profiles which include buoyant stability effects on 
mean velocity, temperature and humidity to a "force-restore” formulation 
for the evolution of surface soil temperature to yield the local values 
of shear stress, heat flux and surface temperature. A self-contained formu- 
lation is presented including parameterizations for solar and infared radiant 
fluxes at the surface. 

In addition to reference- level mean flow properties parameters needed 
to impliment the scheme are the thermal heat capacity of the soil per unit 
surface area, surface aerodynamic roughness, latitude, solar declination, 
surface albedo, surface emissivity and atmospheric transmissivity to solar 
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radiation. 


Sample calculations are presented for a case with constant atmospheric 
forcing at the reference level and for a variable atmospheric forcing case 
at conditions corresponding to Kahle's (1977) measurements of windspeed and 
air temperature and radiometer soil surface temperature measurements under 
dry vegetatlvely sparce conditions in the Mohave desert In California, USA. 

The latter case recovered the observed ground temperature variation over a 
diurnal cycle reasonably well for the parameters used, and displayed a variety 
of buoyant stratification conditions which can occur in atmospheric sur- 
face layers Including convectively unstable, stable and stable-decoupled 
zones. 


1, Introduction 

For a number of applications in micrometeorology and air pollution 
transport analysis it is desired to know the turbulent shear stress x and 
the turbulent heat flux F in the so-called constant flux atmospheric surface 
layer immediately above the earth's surface when direct measurements of the 
relevant turbulent fluctuation moments are unavailable. What is often available 
is information on mean horizontal velocity, temperature and humidity at some 
reference level sj in the surface layer where a time series has been measured by 
appropriate instrumentation. Generally the reference-level mean velocity and 
temperature uj(t) and Tj(t) are time-dependent over diurnal cycles, where the 
overbar average is understood in the usual sense to denote averaging over 
short period turbulent fluctuations. The problem to be considered here is the 
generation of consistent values of i(t) and y{t) from the reference-level mean 
flow using a physical model for the structure of the atmospheric surface layer 
and the underlying layer of ground. 
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An Important consideration In this context is that the aerodynamic 
drag coefficient Cp = t/(p m^ 2 ) and the heat transfer coefficient c H = 
-F/(pCpU 3 (r 3 -T a )) , where p is the air density, C p is the constant-pressure 
specific heat of air and T a is the surface temperature (overbars are dropped 

a 

here and henceforth on reference-level properties), are not constant but 
can depend strongly on the buoyant stratification of the surface layer. 
Generally both C D and C H increase for unstable (heated from below) conditions 
and decrease for stable (cooled from below) conditions. Under extremely stable 
conditions, which are nevertheless encountered in practice, the surface 
layer turbulence can be extinguished entirely, leading to decoupling of the 
surface layer from a possibly turbulent zone some distance above it. Under 
these conditions the drag and heat transport drop substantially to zero, 
laminar transport being negligible in the context of air-surface interactions. 

The structure of the surface layer has been extensively treated in 
recent years in the context of Monin-Obukhov similarity theory in which the 
horizontal velocity and temperature profiles u(z) and ¥{z) are uniquely defined 
by the aerodynamic roughness length z Q , the friction velocity u* = (x/p) 1 ^ 2 , 
the buoyant temperature scale T* = -F/(pC p w*) and the surface temperature 
T g r. T{s q ) (Monin and Yaglom, 1971, p.430). In particular, the Monin-Obukhov 
length scale L = T 8 u* 2 /(KgT *) , where k is von Karman's constant and g is 
the gravitational acceleration, is defined with. dimensionless mean velocity and 
temperature gradients in the surface layer given by the "universal" 
functions $ m (s/L) n (kz/u*)%u/%z and $ H {z/l) ~ [Ks/T*)dT/Zz. In practice, the 
<f>^ and functions are found experimentally and the velocity and temperature 
profiles determined by integration of the gradient functions from z = z 0 
to some arbitrary point in the surface layer. In developing the details of 
the present scheme, the functional forms of $ m (z/L) and 4> fl (s/L ) resulting 
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from an experiment performed In wheat fields of Kansas, USA, by Buslnger et al. 

(1971) are used. It Is recognized that other workers have measured different 
functional forms, for example Hicks (1976) and Slsterson and Frenzen (1978) 
whose results differ primarily In the stable case. In addition HI cl:s (1976) 
has questioned the Inequality between the neutral values of and $ g In the 
Buslnger et al functions which he finds difficult to correlate with the 
physical processes Involved. 

We hasten to point out at the outset that other forms of the Monin- 

Obukhov functions can be used In the framework of the scheme developed here 

* 

as experimental discrepancies are ultimately resolved. Also, we should, 
strictly speaking be working with potential, rather than physical, tempera- 
ture, but for surface- layer reference levels of the order of ten meters or 
less, the results should be insignificantly affected. Humidity, on the 
other hand, can significantly affect the computed fluxes, and its effects are 
Included in the subsequent formulation; it is omitted here only to more rapidly 
focus on the central problem. Namely, given the aerodynamic roughness z c , 
the reference-level velocity and temperature uj and Tj and the surface 
temperature the flux parameters u* and T* are uniquely defined by the 
Monln-Obukhov similarity profiles; however, while s Q can often be characterized 
by the known properties of the terrain in question, and Uj and Tj are known, 
the surface temperature T 8 is almost never measured on a routine basis. And even 
if the soil surface temperature immediatly below the reference-level instruments 
were known, its significance might be unclear since it could well represent 
only the characteristics of the immediate surface type, rather than the regional 
average of interest for flux estimations. 

A more useful approach in practical situations is to characterize the 
thermal and radiant properties of the soil surface layer and to compute the 


5 


soil temperature simultaneously with the Monln-Obukhov parameters at 
each timestep by solving the soil heat conduction equation under the 
Influence of solar and Infared radiation at the surface as well as the 
turbulent sensible and latent heat flux associated with the atmospheric 
surface layer. A number of approaches to the calculation of ground temper- 
ature exist In the literature, primarily in relation to parameterization 
of fluxes in general circulation models, but In all of this published work 
the atmospheric fluxes are represented by constant values of the drag 
and heat transfer coefficients. Perhaps the most accurate treatments of the 
ground temperature unsteady heat conduction equation are the f ;nite-difference 
solutions of Benoit (1976) and Kahle (1977). In the interest of conserving 
computer time a number of integral approaches leading to ordinary differential 
equation models for surface soil temperature have been proposed In recent 
years. A method developed independently by Arakawa (1972) and by the British 
Meteorological Service (Corby et al ., 1972; Rountree, 1975) utilizes a rate 
equation for a grounc "slab" temperature T 8 dependent upon forcing by the 
sum of atmospheric and radiant energy fluxes; however, as noted by Bhumralkar 
(1975) this method omits the influence of soil heat flux on the underside of 
the integral slab. Moreover, the temperature being predicted is not truly the 
surface temperature consistent with the similarity atmospheric surface layer, 
but some depth-average of the temperature in the thermally active layer of 
soil beneath the surface. 

The scheme to be developed here is based on the so-called "force-restore" 
ordinary differential rate equation independently proposed by Bhumrahlker (1975) 
and Blackadar (1976) and recently evaluated in comparison with a 12-layer 
finite-difference solution to the partial differential heat conduction equation 
by Deardorff (1978). The force-restore formulation has the advantage of predict- 
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Ing the soil temperature T at the surface, rather than a depth-aveaged 
value, and contains a mechanism by which a deeper soil layer can Influence 
the surface temperature. Deardorff (1978) found the force-restore formulation 
was computationaly more efficient than the use of multiple soil layers 
and superior to five other approximate methods In current use when diurnal 
forcing was present. In a related study, Deardorff (1977) has proposed 
an analogous force-restore formulation for the soil-surface moisture fraction 
as a parameterization of ground-surface moisture content for use in atmospheric 
prediction models. 

In what follows, the relevant aspects of Monln-Obukhov similarity 
theory are developed, the force-restore equation for soil surface temperature 
Is derived, expressions are given for the solar and Infared radiation forcing 
terms, and the coupled system computer code Is described. Subsequently, 
computational results are presented for both constant and variable atmospheric 
forcing and comparisons made with observational data. 

2. Monin-Obukhov Profiles in the 
Atmospheric Surface Layer ' 

The planetary boundary layer (PBL) is the region of the lower atmosphere 
where flow is turbulent on a microscale by virtue of shear- and convection-induced 
turbulent fluctuations associated with the underlying solid or water surface. 

In the lower PBL, i.e. the so-called surface layer where the turbulent shear 
stress, turbulent sensible heat flux and turbulent humidity flux may be 
treated as constant, the relevant flow variables are velocity (u, v, w) = 

(u + m ', o ', !.»*), temperature T = T + T' and specific humidity (7 = q + q ' , 
where overbars denote Reynolds aveages and "primes" denote fluctuations. We 
assume for the present problem that the known reference- level properties Uj, 

Tj and (7^ represent u t T and </ at some height Zj within the surface layer. 
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What we would like to know however are the Reynolds stress t * -puV , 
sensible heat flux F * pCpW'I” and evaporative flux E = p i/q r which may be 
considered constant In the surface layer. The mean air density p may also be 
treated as constant at Its surface value ( e 1.23 kg/m 3 ), but fluctuations In 
density are Important In buoyancy-generated (convective) turbulence. It follows 
from the equation of state p * p/[R(l + 0»61^)r that the effects of water 
vapor variations on density can be handler by working with the virtual temp- 
erature T v = (1 + 0.61<7)3’. Accordingly, the Bousslnesq approximation (which 
neglects pressure fluctuations) for air with some water vapor relates density 

fluctuations to fluctuations In virtual temperature, p'/p a -Ty'/iy Ordinarily 

* 

we can take T - T v except when differences or fluctuations are involved. An 
Important quantity In buoyantly-interactive turbulent flows is the buoyancy 
flux Fj, = -CpF’pV “ pC’pW'i’y' = F + 0.61C p TE t where E t F and F b are all pos- 
itive upwards, and buoyantly unstable, neutral and stable conditions corres- 
pond to F b negative, zero and positive, respectively. 

To relate the mean flow measurements in the surface layer to the surface 
layer fluxes we shall make use of the similarity theory orginally proposed 
by Monin and Obukhov and subsequently developed by many others. 

It is convenient here to introduce turbulent velocity, thermal and humidity 
scales related to the surface fluxes as follows, 

u* 5 (t/p) i/2 , T* - -F/(pC p u*) t q* = - E/(pu *). (2) 

In addition, the buoyant fluctuation thermal scale 

V 5 - F b f (* c p u *) = T * + 0.613V 4 (3) 

plays an important role in Monin-Obukhov theory. The fundamental assumption 
based on dimensional analysis arguments is that vertical gradients of mean 


- 8 - 



( 4 ) 


flow properties In the surface layer may be expressed, 

3u u** m (z/L) 3 f T** H (3/L) 3 q q** q {z/L) 

iz K3 3S K3 dS K3 

where the coefficient tc known as von Kantian' s constant has a numerical value 
of k “ 0.35 (Buslnger ei at . « 1971), 

I = T g u* 2 / (<gT v *) (5) 

Is a lengthscale characterizing buoyancy effects on turbulence known as the 
Monln-Obukhov length, and ^ and ^ are "universal" functions of the 
surface layer stability variable c = z/L. In general L Is negative for unstable 
stratification, positive for stable stratification and approaches positive dr 
negative infinity under neutral conditions. Businger et at. (1971) have measured 
and curve fit the fuctions <f> m (c) and $#(?) under both unstable and stable 
conditions. They find 

*„(c) - (1 - ,„({) . a 0 (l - 9c)- l/2 

under unstable (c<0) conditions and 

t m U) = 1 + 4.7c . $ H U) = «o + 4 -7c 

under stable (c>0) conditions, where the turbulent Prandtl number a = has 

the numerical value at neutral stability of a Q = 0.74. More generally, of 

course a is itself a function of c* On the other hand Dyer(1967) has shown 

<t> /$„ a 1 under a range of unstable conditions. Accordingly, we shall take 
q U 

(c) - (c) In the following analysis. 

H 

Now, integrating the relations in (4) between the "surface", located at 
aerodynamic roughness height s 0 where u(z 0 ) = 0, T{z 0 ) = T e and q{z Q ) - q B 
by definition, and some arbitrary height z ^ in the surface layer, gives the 
Monin-Obukhov profiles for turbulent-mean velocity, temperature and humidity: 
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u[aj) - (u*/ic)[An(a J /a 0 ) - t u (sj/£)], 
Ti*l) " T e + [^(ao/KJltmtaj/a,,) - ^(aj/i)], 
?(» 3 ) " ? 8 + [<? 4 (a 0 /K)][in(a J /a 0 ) - * 7 (aj/£)] ( 


( 6 ) 

(7) 

(8) 


where the stability-dependent profile functions \p u and 1> T are defined 
formally by the Integrals, 


. , . ,h n-VdJ* , C J 

“/ * 

5 0 5 

♦ Uj/l) i / t: [!:;? Tl V (c)3dc . /'i [ _ ° . 

' 6 ' 

These profile functions may be obtained explicitly, for example, by substi- 
tuting the $ m (c) and expressions of Buslnger et al. given previously 

and integrating. The result (obtained after some algebra) Is summarized In 
Table 1. The- mean profile functions (6)-(8) with the stability-dependent 


Table 1. Stability-dependent functions in Monin-Obukhov profiles. 


stabil ity 
variable 

\ * 

* u (a 2 /L) 

tyfizj/L) 


f T+Cl-15(a 7 /i)]l / ^ 

H 2 ) 


(zj/L)<0 
[unstable ] 

H 2 ) 

f l +[l-9(rj/L)]l/2, 

H 2 ) 


- 2ton- 1 [l-15(a J /l)] l/4 + v/2 


(*JL)> 0 
[stable] 

- 4.7(zj/L) 

| 

- 6.4 {z 2 /l) 
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functions of Table 1 are substantially the same as the surface layer 
profiles used by Deardorff (1972) In his PBL parameterization scheme. 

Now, assuming these stability-dependent profiles describe the vertical 
variations of velocity, temperature and humidity In the surface layer, we may 
apply them to the present problem of finding the surface layer fluxes as 
follows. Note first from (3) and (5) that tfv- Monin-Obukhov similarity 
length scale can be written 

u* 2 T 

L = * t 

K g(T* + o.6i:r e q*) 

where the gravitational acceleration g Is approximately 9.81 m/s 2 at sea 
level. Using equations (6)-(8) to express the turbulent velocity, temperature 
and humidity scales u* t T* and q* t and substituting these into the o^jve 
leads to the following Implicit equation for Monin-Obukhov length: 

m h £ („ • 

^[r r r e +0.61r e (q r q e )] [ln{z 2 /z 0 ) - i> u (z 2 /l)Y 

Accordingly, if the reference-level velocity u 2 s u(z 2 ), temperature T 2 = r(sj), 
and humidity q 2 a q(z 2 ) are known, for example from surface observa- 

tions, the aerodynamic roughness z Q is specified, for example as a function 
of the underlying terrain, and if the "surface" temperature T g - f{z Q ) and 
humidity q g = q(z 0 ) at the ground/air interface are known, for example from 
a simultaneous heat transfer analysis as described later in this report, then 
the Monin-Obukhov length L is gi' an in principle by the solution to (9). For 
the unstable (L< 0) case the ty u (z 2 /L) and ty T {z 2 /L) functions are transcendental 
(Table 1) and an explicit solution for I is not feasible. However, experience 
has shown that a straighforward Newton-Raphson iteration to find the root of 
f(L) * 0 is generally effective in the unstable case provided the initially 


- 11 - 



guessed value of L has the proper (negative) sign. 

Uncer stable (z>0) conditions the profile functions are much simpler, 
♦ U («j/L) * - 4 IzjjL and ■ - 6.4*j/L (Table 1), and equation (9) reduces 
to a quadratic In L, 


aL z + bL + o a 0, 


where a * [£n(3j/a 0 )] 2 , b * in[zj/z 0 ) (9.4 - a Q u 1 2 T g /{giT 1 ’-T fj +0.b]T 8 (q 1 -q 3 )l )) 

and a - 4.7aj(4.72j - u 2 2 T B /{g[T 1 -T e +Q.b'\T e {q 1 -q 8 )']}) . Notice here that we 

can Identify the quantity T-T +0.613* (<? 7 -q 0 ) with the difference In virtual 

temperature between the reference height and the surface T^-T vg . This is 

to be expected since virtual temperature differences are related to buoyancy 

differences In humid air. The solution to the quadratic takes the u:ual form, 

- b + Jb 2 - 4 ac 
L = , 

2a 

for the stable (l>0) case, although the nonphysical negative root has been 
discarded. For a real positive root to exist we need 4 ac < 0. Since the coeffi 
cient a Is always positive [being the square of £«( 2 j/ 2 0 )], this means o 
must be negative , that is:4.72 7 - «. 2 T /[g{T -T )] < 0. This condition may 
be written in the form 


W 


< 0.21 


which is basically a condition on a finite-difference form of the Richardson 
number Ri = (gZT v /Zz)/[T 8 (*u/2z) 2 '] corresponding physically to the requirement 
that the Richardson number must be below some critical value Ri < Ri „ = 0.21 
for turbulent flow to exist under stably stratified conditions. In practice 
this means surface and reference-level properties must correspond to Richard- 
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son numbers below 0.21 for a real positive Monin-Obukhov length to exist and 
be calculable for the stable case. Generally speaking* a Richardson number 
above 0.21 computed in this way would indicate a reference level above the 
presumably shallow surface layer, and an associated decoupling of the surface 
turbulent zone from mixing zones above. 

To summarize, given the aerodynamic roughness z Qt the reference-level 
measurements u lt T 2 and ^ and the surface temperature and humidity T and 
q 8 , we can find the Monin-Obukhov length L from equation (9). Generally, we 
would expect L to vary markedly with time-of-day, with positive stable values 
characterizing the nocturnal ground-based inversion, a transition though an 
infinite neutral condition near sunrise, and a subsequent evolution of negative 
unstable values during the day with a negative minimum in mid-afternoon 
corresponding to peak ground temperatures, followed by a transition back Lu pos- 
itive stable values near sunset. Knowing £, the frictional velocity, temperature, 
and humidity scales ca.i be found from 

w* = * UjUtnizj/sJ - (10) 

T* - K{T 1 -T 8 )l[<x 0 [Zn{zi/z o ) - ^(z^/i)]) , (11) 

q* = <{q 1 -q 8 )/{(i 0 iln{z 1 /z o ) - * y ( z^L )]) = T* ; % q j-q g ) / (T j-T g ) , (12) 

which follow immediately from equations (6)-(8). 

3. Surface Soil Temperature from 
Force-Restore Rate Equation 

Since T e and q s are not specified by the surface observations, it is 
necessary to introduce some additional considerations to evaluate them* A s 
discussed in the Introduction, our approach to computing T 8 and q 8 at each 
timestep is to use- the force-restore formulation for the soil surface temp- 
erature (Bhumralkar, 1975; Blackadar, 1976; Deardorff, 1978) as developed below. 
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Strictly speaking T Is the air temperature at the a level of the 

3 O 

atmospheric surface layer. We shall assume however* that air at this level 
Is In thermal equilibrium with the underlying solid or liquid surface* neg- 
lectlng any possible effects associated with Intermediate layers. Accordingly, 
we can identify T with the ground (or soil) temperature at the surface. 

o 

To a good approximation this temperature Is governed by the one-dimensional 
heat conduction equation for a semi-infinite slab heated or cooled at the 
surface by radiation and the turbulent sensible and latent heat fluxes. By 
the Fourier heat conduction law the rate of heat flow vertically through the 
soil at depth a below the surface Is proportional to the temperature gradient. 


G(s,t) * - xar/33. 


where X is the coefficient of thermal conductivity; a function in general 
of soil type and water content. A typical value Is X = 2. 5x10" s cal/ # K*cm*s 
(- 1.05 J/°K*m*s) for soil, although Fig. 37 of Sellers (1965) indicates 
variations of a factor of two or more from this value are certainly possible 
depending on soil type and water content. 

Since the soil surface temperature varies largely in response to 
diurnal cycles of solar radiation we need to consider the unsteady form of 
the soil heat conduction equation, 

3Q 3 2 T (13) 

— = = X — , 

3fc 3 xj 3a 2 

where p is the soil density, a is the soil specific heat per unit mass. Clearly 
then, the product p<? represents the soil's specific heat per unit volume. The 
quantity K = \/(p o) may be defined as the thermal diffusivity. Here a repre- 
sentative value for soil is K - 5x10 3 cm*/s (5x10 7 m 2 /s), although again 
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Sellers (1965, Fig, 39) Indicates factor-of-two or more 

variations in K are possible depending on soil composition and water content. 
Accordingly, a typical specific heat per unit volume of soil is 
p o ■ \/K • 0.5 ca1/°K*cm 3 («*2.1*10 6 J/°K*m 3 ). 

Noting that the vertical coordinate In equation (13) is positive down- 
ward, the relevant boundary conditions on (13) for the present problem are 

G(0,t) - -A(ar/3a) a - + 5 - /? - F - lE t (14a) 

<?(«.) - 0, (14b) 

\ 

where s Is the net flux of solar radiation absorbed by the ground surface, 

J? is the net longwave flux radiated up into the atmosphere, F Is the sensible 
heat flux and Lff is the latent heat flux leaving the surface by turbulent 
transport vertically through the atmospheric surface layer. In the general 
case the net heat flux to the ground at the surface &'(0,i) is too complicated 
a function of the time t to admit a simple, analytic solution to (13) with 
these boundary conditions. We can recognize however that surface temperature 
variations over timescales of , say, a few days are driven largely by the diurnal 
solar radiation cycle of frequency ft 3 2* radians/day - 7.27*10“ 5 s -1 . It Is 
Instructive therefore to consider the response of (13) to a periodic boundary 
condition of the form 


r{Q.t) « <'<v + AT^itiM(at). (15) 

Assuming a semi-infinite solid below, Carslaw and Jaeger (1959) show that the 
solution to (13) with this boundary condition has the fonn 

VOM) * <V + A ry a/ti iu»(M - */</), (16) 

where 
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d = (ZK/n) l/2 = (2x/ (pen)) 1/2 (17) 

has the significance of a diurnal skin depth for the penetration of a 

thermal wave of period n applied at the surface. Thus, the typical thermal 

diffuslvlty for soil of K - 5*10” 7 m 2 /s quoted earlier corresponds to a 

characteristic diurnal skin depth of some d = 0.12 m (12 cm). 

The force-restore approximation Is based on the fact that while the 

surface forcing terms of (14a) are not precisely periodic of a sinusoidal 

form, the influence of the solar forcing of period n establishes the penetration 

of diurnal waves in the manner of the exact solution, altough the actual 

boundary condition on the heat flux of (14a) must.be satisfied to transfer 

% 

the proper amount of heat to and from the soil. Taking the partial derivative 
of (16) with respect to z and multiplying the result by -d gives the relation, 

-d(%T/Zz) = LT 0 e -3 ^(c?os(nt - z/d) + sin(ut - z/d)) . 

Moreover, taking the partial t derivative of (16), multiplying the result 

* 

by 1/n, adding T from equation (16) and subtracting <T> from both sides of 

o 

the equation leads to the relation, 

1 87 / , 

-• — + T - <T S > - AT ' O e~ 3/ci (cos(nt - z/d) + cin(Qt - z/d)). 
n 3 1 v 

Recognizing that the right-hand-sides of the foregoing two equations are 

identical, and that the heat flux through the soil at any depth is G = -X3T/33, 

we may eliminate the group of terms on the right-hand-sides to get 

1 3 T 3 T d( 3T] d 

+ T - <T S > = -d - -A— = -G. (18) 

SI 3/ 3a A l 3sj A 

The force-restore ordinary differential equation for surface soil temp- 
erature is now obtained by evaluating (18) at the surface (a = 0, T = T ) 
using the actual (nonsinusoidal , in general) boundary condition (7(0, t) 
from (14a): 
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( 19 ) 


dT e 2 

— - — (S - K - F - L •&•) - n(r - <r„>), 

dt c e 

where the heat capacity per unit surface area appearing above 


c 8 = P = 2x/(n d) (20) 

corresponds to diurnal forcing cycles. 

In the periodic solution <T> as defined in (15) is the mean soil 
surface temperature; in the force-restore equation (19) it was suggested 
by Blackadar (1976) that for short range projections <T S > be treated as 
a constant whose value is estimated from the mean air temperature for 
the prior 24 hours. For forecasts of three days or more Deardorff (1978) 
suggests the variation of <T 0 > be computed from 


d<T > 1 

= (S - R - F - L ♦£), (21) 

dt <c t > 

where =• (365) l/2 c’ g is the heat capacity per unit surface area for 
the annual thermal wave. 

In the sample calculations presented later, we are interested in 
diurnal cycles, so <T > is treated as a constant arid c, has the diurnal 
value. For the prior typical soil values of pe = 2.1xl0 6 and d - 

0.12 m the heat capacity per unit surface area is C 8 c 2.5x10 s J/°K-m 2 . 

Notice that in the force-restore formulation (as in the integral formulations 
of Arakawa (1972), etc.) the various thermal properties of soil are sub- 
sumed into the one parameter <’ J} . Notice that the last term in (19) acts 
to restore T s exponentially to some mean (or deep) soil temperature if the 
surface forcing terms are removed. To use this differential equation in 
the present context we need to express the fluxes H, F and K in terms 
of the local values of r <: , uj, T } , and parameters characterizing the 
parti cl ar site. 
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Consider first that the sensible heat flux Is expressible from (2), 


(10) and (11) as 


F = -p C p u*T* = , (22) 

where 

C H 3 C l = K 2 /(a 0 [tn(3 J /3 <7 )-^ u (3 J /ii)][tn(3 J /3 o )-i|/ 2 ,(a J /ii)]) 

Is the surface layer's heat transfer coefficient referred to level Zj. Since 
this flux Is positive upwards we generally have f>0 when Thus, know- 

ing T 1% T g and q e we would first find L from (9); the corresponding 
sensible heat flux follows immediately from (22). 

Analogously, the latent heat flux carried by water evaporation or 
condensation at the surface is expressible from (2), (10), (11) and (12) as 

IE = -pLwV = ~C H 'plu 1 (q 1 -q c ). (23) 

In principle this can be handled similarly to the sensible heat flux although 
a major problem still exists insofar as we have not yet specified how to 
find the surface humidity q g . A number of approaches exist, all of which 
require additional consideration of evaporation and evapotranspi ration in 
the soil-vegatative component of the hydrological cycle. However, a fundamental 
idea in all of these methods is that the actual evaporative rate cannot 
exceed the potential evaporation rate which would obtain if the humidity at 
the surface were saturated at the value corresponding to T s , 

E o = ~ C H’ pu l^l^ 8a t {T ^' (24) 

where q gat .{T ) is calculable from the Clausius Clapeyron equation. 
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It Is useful sometimes to represent the actual evaporation 
rate In terms of an actual -to-potential evaporation ratio 0 = E / E 0 As' 
discussed by Sellers (1965), It Is helpful also to distinguish two stages 
of evapororatlon from soil and vegetation which depend on the volume fraction 
of soil moisture In the active soil layer w, usually expressed In millimeters 
of water. In the first stage, when the soil moisture content is greater than 
some critical value W evapotranspiration proceeds at about the potential 
rate E Qt and depends mainly on external meteorological factors (0 * 1 when 
!/:> Wj,). In the second stage, when the soil moisture content is less than the 
critical value, the rate of evapotranspiration depends on the soil moisture 
content, with the relationship often assumed to be linear (0 = w/w^). 

Clearly, when the soil is below the saturation value w k , it is necessary to 
either measure or model the evolution of W to predict evaporation. Deardorff 
(1977), for example, has proposed an extension of the force-restore approach 
as a prognostic equation for w driven by the difference between evaporation 
and precipitation rates and a restore term proportional to the difference 
between the local W and a long-term average <W 8 >. Thus an effective 0 can 
be computed simiiltanuously, or simply prescribed. Note that equating the actual 
evaporation rate -Cjf'Wjiqj-q ) to $E o gives a relation for the surface 
humidity in terms of 0 , q 2 and T g , 

q s = (1-ekj + (<7 y ± w)- 

In turn, this relation can be used to specify the surface specific humidity 
appearing in equations (9), (12) and (23). 

Turning now to the solar radiation term 5, we recognize first that 
diurnal variations in heating are associated with the daily variation in the 
solar zenith angle 7. (the angle between a line pointing toward the sun and 
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a vertical line normal to the earth's surface at the latitude of Interest 
— e.g., when the sun Is directly overhead Z « 0°). It can be shown that 
the cosine of the zenith angle Is expressible as 

cosZ(t) = 8tn$8in6 + cosQcoe 6ao8 [sit ) , (26) 

where $ Is the latitude, 6 Is the solar declination angle and fit Is the 
hour angle relative to solar noon (i.e, nt = 0° @ 1200 hr and Increases by 
15° every hour so that, for example, fit * -90° @ 0600 hr). The solar declina- 
tion'. varies slowly relative to the hour angle In a sinusiodal fashion with 
a period of one year and a maximum amplitude at th.e summer solstice (June 
21) of 23°27' and minimum at the winter solstice of -23°27'. Values for 
each hour and day of the year may be obtained from The Nautical Almanac 
published by the U.S. Government Printing Office although for the present 
analysis we may justifyably take <5 constant over a diurnal cycle. It 
should be clear, however, that the zenith angle only has significance 
during daylight hours when z lies between 0° and +90°. The hour angles 
(or times of day) corresponding to local sunrise or sunset are therefore 
found by setting cosz = 0 in (26) and solving for the hour angle, 

cos (fit) = - tan$tan&. 

The two possible values of fit between -90° and +90° correspond to local 
sunrise and sunset, repectively. 

Under cloud-free conditions, the direct solar radiation absorbed by 
the ground may be written, 

S(t) = S 0 co8Z ( 1 - A)x 8ecZ , 

where S 0 is the solar constant, A is the surface albedo and t is the atmospheric 
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transmission coefficient. The commonly accepted value of the solar constant 
(the frequency-integrated solar radiation flux per unit area falling on 
a plane perpendicular to the sun's rays at the top of the atmosphere) Is 
S Q « 2.0 lymln” 1 (2.0 cal •m1n” 1 *cm“ 2 « 1400 W*m -2 ) which we adopt here 
as well, for numerical calculations. Sellers (1965, p. 21) has tabulated a 
range of albedos, or reflectivities, of various surfaces in the shortwave 
portion of the electromagnetic spectrum (wavelengths less than 4.0 urn, the 
region relevant to reflection of solar radiation). By definition, these 
albedos are between 0 and 1 and increase with increasing reflectivity; thus 
we have A in the range .05-. 15 for coniferous forests, .1-.2 for deciduous 
forests and green meadows, .15-. 20 for tundra, chaparral and wet-season 
savanna and .25-. 30 for dry-season savanna and desert. In the present study we 
adopt a nominal value of A - 0.12 for purposes of Initial calculations. 

For the transmission coefficient t, values In the range of 0.75-0.90 are 
typical of the fraction of solar radiation penetrating to the surface under 
cloud-free overhead sun conditions. It is worth noting here again that the 
time dependence of s(t) in the cloud-free case is dominated by the Sit 
diurnal periodicity of the zenith angle consistent with the assumption 
made earlier of an effective diurnal penetration depth for the solar heating 
wave. 

A semi -empirical correction for solar radiation can be developed for 
overcast skies in terms of the degree-of-cloudiness parameter n often available 
with the surface observations. We note first that the clovd-free term 
S* = cop?. y BC(%7 ‘ represents the direct component of radiation from the sun 
through the atmosphere incident on the surface. Some fraction of this is 
diffusely scattered by the atmosphere and reaches the surface as well, although 
this is usually neglected or treated implicity in the cloud-free case. The 
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fraction of diffuse to direct solar radiation e B S^^fs* is solar zenith 
angle dependent varying from 5% for an overhead sun to 15% for the sun at 
the horizon. Assuming a simple linear variation with oosZ (Kahle, 1977) gives 

e(z) *= 0.05 + 0.10(1 - oobZ ). (27) 

Now, under cloudy sky conditions, the direct component of radiation is 
reduced to S*(l - rt) as discussed by Kondratyev (1969, p.312), while the 
diffusive component becomes more important since it now Includes the 
diffusive scattering by clouds. Kondratyev (1969, p. 399) writes a parameter- 
ized form of the diffuse radiation flux under clpudy conditions as 
S*fe(l - n) + tf«(i + e)], where AT is an empirical latitude- dependent param- 
eter representing the solar radiation transmitted by diffuse radiation 
through clouds. Based on Table 8.5 of Kondratyev (1969, p. 468), however, 
the variation of K with latitude is relatively weak in the middle latitudes, 
with numerical values In the range of 0.32 to 0.36 from t = 0° down to 
$ B 55°. For our calculations, the constant value K = 6.34 Is adopted for the 
ratio of diffuse radiation transmitted through clouds to the direct component. 
Accordingly, the total solar flux incident vertically on the ground In the 
presence of clouds may be written s*[(l - n) + e(l - n) + *n(l + e)]. Re- 
arranging, and using the numerical value of K, we may wri „e the solar flux 
actually absorbed by the ground In the form, 

s(t.n) = S 0 cooZx scoZ (1 + e(z)) (1 - >1 ) ( 1 - 0.66*), (28) 

where e(z) Is given by (27) and z(t) by (26). 

The remaining term on the right-hand-side of (19) to be parameterized 
is the terrestrial longwave (infrared) radiation R upwelling from the surface. 
Observations indicate that this flux is correlated under cloud-free conditions 
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with the ground surface temperature T g and the water vapor pressure a 
a few meters above the surface, say B e(*j), by expressions of the . 
form R ■ or e 4 /(ej)» where a Is the Stefan-Bol tzmann constant * 8.14 x 10" 11 
lym1n _1 *°K“ 4 (8.14 x 10" 11 cal^crn^-mln" 1 *^" 4 - 5.86 x 10' 8 W*nr 2 *°K" 4 ) and 
the function f(e 2 ) Is determined empirically. For the present model, we 
use the formulation for net upwel ling surface Infrared radiation In the 
cloud-free case proposed by Brunt and quoted by Sellers (1965, p. 53) In the 
form 


R ■ eo7 8 4 ( 1 - a - b/e^) % 

* 

% 

where e, here. Is the surface eneilslvlty and a and b are empirical coefficients 
Ordinarily the vapor pressure In this expression Is expressed In mm Hg. This 
can be computed readily from the reference-level specific humidity q 2 (kg/kg) 
and the station pressure p (In. Hg) by 

ej(mm Hg) = 1.61<?j x p(1n. Hg) x 25.4 mm/in. (29) 

Sellers (1965) in his Table 7 quotes surface emmi si vi ti es for varl 
ous natural and vegetative surfaces in the range of 0.9 to 
0^95. For. the coefficients In Brunt's formula we take Budyko's 
(1956) values, a - 0.61 and b B 0.050 (mm Hg)" 1 / 2 . These give results sim- 
ilar to the values In Kondratyev's (1969) book and are close to the median 
of twenty-two evaluations quoted by Sellers (1965). In reality, however, the 
sky generally contains some cloudiness, and the Infrared flux term should 
be corrected approximately to account for this effect. Kondratyev (1969, pp. 
575-576) suggests an empirical correction of the form R = /?*(1 - cn), where 
R* Is the cloud-free value of infrared radiation, a => 0.76 is an empirical 
coefficient, and n Is degree-of-cloudiness. Notice that clouds have a blanket- 
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Ing or Insulating effect on ground temperature Insofar as they reduce 
the radiant heat loss from the surface through back radiation. Combining 
the relations and numerical factors derived thus far gives the following 
parameterization for the terrestrial Infrared term In (15), 

R(T a ,q v n) « 0.92or s 4 (0.39 - 0.050^7)(1 - 0.76n), (30) 

where Is related to qj by (29) . 

4. Numerical Model and 
Sample Results \ 

% 

The foregoing scheme for generating surface layer turbulent fluxes 
from mean flow observations has been impllmented In a FORTRAN computer 
code (TURBFLUX) which has been run thus far on the CDC 6600 machines at 
Brookhaven National Laboratory and the Courant Institute of Mathematical 
Sciences at New York University. The logic of the calculation is reviewed 
below. 

The basic Input data are ti.e reference- level wind speed u^(t), temper- 
ature relative humidity rj[t) and fractional sky cover n(t) in dig- 

ital form over the period of interest, the reference-level height sj, 
aerodynamic roughness s 0 , specific heat per unit surface area c , initial 
surface temperature r s (0), mean soil temperature <T e > , latitude <{>, solar 
declination 6, surface albedo A , surface emissivity e, ground wetness par- 
ameter p and integration timestep Af. At the initial time, and at all subsequent 
timesteps after T is computed, the following routines are execute 

(1) Compute <7j and from r^, T lt and equation (25). 

(2) Compute T v j and T vg and evaluate Richardson number, Ri j - psj(r y ^-r t) J/( 

(3) If decoupled (0 ,2\<Rij) set u* = T* = q* - 0, If stable (0<A't j<0.21) find 
l from quadratic equation. If unstable [Rtj< 0) find l from Newton-Raphson 
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Iteration solution to equation (9). 

(4) From uj, Tj* qj* T a and q g find u* t T* and q * for stable, neutral of 

stable cases from equations (10), (11) and (12) and Table 1 functions. 

(5) Compute F from u*, T* and equation (22). 

(6) Compute L«E from «*, q* and equation (23). 

(7) Compute S from t, n and equation (28). 

(8) Compute R from T 8* <1 2 and n from equation (30). 

These operations completely define the right-hand-side of the force-restore 

equation (19) at each tlmestep. To improve the accuracy cf the integration 

% 

a semi-implicit technique Is used to evaluate T g {t) numerically from (19) 

In the TURBFLUX code. For the calculations discussed next an integration time 

of Ai = 10 min (600 s) was used with outputs printed every hour. 

The first case studied was for constant atmospheric forcing, dry soil, 

mld-lattltude equinox conditions for the soil and atmospheric parameters given 

« 

In the caption of Figure 1. Here, the reference-level windspeed and temperature 
were held constant at = 4 m/s and Tj - 280 °K over a diurnal cycle. As 
shown In Figure 1(a) the friction velocity varies slightly about Its neutral 
value, being somewhat higher during the unstable daytime phase and lower 
at night. The transition from unstable to stable surface layer flow is marked 
by the sign change In the buoyant temperature scale T* and occurs slightly 
after sunset, with another transition back to unstable turbulent flow the 
following morning. This relatively smooth variation in the buoyant stability 
of the surface layer can also be traced in the variation of the reciprocal 
Monin-Obukhov length shown in Figure 1(b). Also shown is the computed variation 
of soil surface temperature. Notice the rapid drop in the afternoon as the 
sclar radiation diminishes, followed by a somewhat slower radiational cooling 
at night with a subsequent build-up the following day. 
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FIGURE 1 Typical diurnal cycle of atmospheric surface laypr with constant atmospheric 
forcing (0 2 ? = '■ 3 m; u; = 4 m*s _1 , Tj = 280.0 °k) and dry soil (6 = 0) 
computed frofii preM-.t model for the soil and atmospheric parameters: 

C g ~ 1 , 6x 10 5 J*nf , s 0 = 4*10"** m, <T,,> = 282.0 °K, $ - 45 deg, 6 - 
0 deg, A c 0.25, e s 0.90 and t = 0.85. (a) Friction velocity u* and 
buoyant stability temperature scale T* and (b) soil surface temperature 
T 8 and reciprocal Monin-Obukhov length ML. 
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In practice, the atmospheric reference- level properties are likely 
to vary significantly over a diurnal cycle. Accordingly, we chose as our 
Initial test of variable atmopsherlc forcing to model the observational 
conditions of Kahle (1977) who measured both reference-level wlndspeed 
and temperature and soil surface temperature (by radiometer) at sites in 
the Plsgah Crater-Lavlc Lake region of the Mohave desert in California. 

This is an arid, vegetatlvely sparce region whose surface 
Includes both basalt and clay playa zones. Soil moisture measurements 
by Kahle (1977) indicated that latent heat transfer could be ignored in 
the surface energy balance. 

The variable reference-level winds and temperatures used in our simu- 
lation of this case are shown In Figure 2, as are the other input parameters 
which are given in the caption. Shown in Figure 2(b) is the model -computed 
surface temperature compared with the radiometer data, where the error bars 
on the observations correspond to a range of values for the region. The 
parameter values of latitude, solar declination, albedo, surface emissivlty 
and transmissivity used in the model are those given by Kahle (1977) and 
the su. ace roughness was chosen to recover the heat transf •' coefficient 
used in Kahle's (1977) model under neutral conditions. The value of C 8 was 
adjusted however to obtain the solid curve, a reasonable procedure in view of 
the semi -empirical nature of this parameter. The complex variations predicted 
for the surface layer with this forcing are, however, better revealed in 
Figure 3. The variability of u* and" 1 * shown in Figure 3(a) are for this 
case the result of variable atmospheric forcing at the reference-level as 
well as variations in stability associated with surface temperature changes. 
In addition, the variation of Richardson number an Monin-Obukhov ler for 
this case are considerably more complex than for the constant-forcing case. 
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FIGURE 2 Diurnal cycle of atmospheric surface layer with variable atmospheric 
forcing (@ z-\ - 1.5 nr, u ^ = uAt ), Tj = TjU)) and dry soil (B * 0) 
in the Pisgan Crater-Lavic Lake region of the Mohave desert in California 
including model results for the soil and atmospheric parameters: 

C = 2.5x10 s 0*m _2 *oK~ 1 , = 3*10 -4 m, <T S > = 301.16 °K (28 °C), * = 

34.65 deg, 6 = 3.1 deg, A = 0.44, r. = 1.00 and x = 0.80. (a) Measured 
windspeed at reference level and (b) Measured air temperature at reference 
level and comparison of measured and computed surface soil temperatures. 




0.2 r- 


1.0 


\ 



FlU'JRt 3 Diurnal cycle of atmospheric surface layer properties computed with 
present model for variable atmospheric forcing and parameters of 
FIGURE 2. (a) Friction velocity u * and buoyant stability temperature 
scale T* and (b) reciprocal Monin-Obukhov length l/r.-and reference-level 
Richardson number Rij. 
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As Indicated In Figure 3(b) the computed evolution of Richardson number Includes 
an unstable region until approximately sundown near 1800 hr, followed by 
an Increasingly stable zone which becomes decoupled at about 2000 hr and 
reattaches as a stable layer around 2300 hr, followed by a transition back 
to unstable conditions which persist into the next day. Notice that the 
"unusual" transition to unstable flow in the noctural phase corresponds to 
the drop In air temperature below the surface temperature slightly after 
midnight, while the large negative values of Richardson number and recipro- 
cal Monin-Obukhov length In the early morning of 30 March correspond to 
the low windspeeds at these times. Of particular interest is the model's 
ability to predict decoupling of the stable layer, which would not 
be possible with a constant heat transfer coeeficient. 

5. Concluding Remarks 

The scheme documented here for finding the turbulent fluxes at the 
bottom of the planetary boundary layer from measurement of mean flow prop- 
erties at some reference height in the range of 1-10 meters above the surface 
has been tested with reasonable success against surface temperature data over 
a diurnal cycle. The method is based on mating Monin-Obukhov similarity 
theory with the- force restore formulation of the ground temperature equation 
developed originally by Bhumralkar (1975) and Blackadar (1976). Indeed 
the present model parallels parts of Blackadar's (1976) model for the 
nocturnal boundary layer, although his reference level is driven by prognostic 
equations rather than direct observations. In view of the current interest 
in understanding and modeling the decoupling of very stable layers at night, 
it would be interesting to test the model's ability to predict such decoupling 
in a controlled observational situation where decoupling actually is measured. 
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